#!/usr/bin/perl
require "grid_utils.pl";

$lcrootdir = shift;
$lctype = shift; # "MODIS" or "NLCD_INEGI"
$lcid = shift; 
$lcpfx = shift; 
$lc_table = shift; 
$domain = shift;
$landmask = shift;
$startyear = shift;
$endyear = shift;
$coarse_tile_mask_file = shift; # can be null
$latmin_in = shift; # southern boundary of region to process
$latmax_in = shift; # northern boundary of region to process
$lonmin_in = shift; # western boundary of region to process
$lonmax_in = shift; # eastern boundary of region to process
$pix_per_deg_y = shift; # modis pix_per_degree, = 240
$resolution_out = shift; # in degrees
$outdir = shift;
$outpfx = shift;
$force = shift;

if ($lctype eq "MODIS") {
  $lcformat = "modis";
  $lcsubdir = "/PFT";
}
elsif ($lctype eq "NLCD_INEGI") {
  $lcformat = "asc";
  $lcsubdir = "";
}

@varnames = ("LAI","NDVI","albedo");

if ($coarse_tile_mask_file ne "null") {

  # Read coarse tile mask
  open(FILE,$coarse_tile_mask_file) or die "$0: ERROR: cannot open $coarse_tile_mask_file for reading\n";
  $row = 0;
  foreach (<FILE>) {
    chomp;
    @fields = split /\s+/;
    if ($fields[0] =~ /ncols/) {
      $ncols = $fields[1];
    }
    elsif ($fields[0] =~ /nrows/) {
      $nrows = $fields[1];
    }
    elsif ($fields[0] =~ /xllcorner/) {
      $xllcorner = $fields[1];
    }
    elsif ($fields[0] =~ /yllcorner/) {
      $yllcorner = $fields[1];
    }
    elsif ($fields[0] =~ /cellsize/) {
      $cellsize = $fields[1];
    }
    elsif ($fields[0] =~ /nodata/) {
      $nodata = $fields[1];
    }
    else {
      for ($col=0; $col<$ncols; $col++) {
        $coarse_tile_mask[$nrows-1-$row][$col] = $fields[$col];
      }
      $row++;
    }
  }
  close(FILE);

}
else {

  $nrows = 1;
  $ncols = 1;
}

# get lists of MODIS tiles to download and process
if ($coarse_tile_mask_file ne "null") {
  # define hvs for each 10x10 tile
  for ($i=0; $i<$nrows; $i++) {
    $lat = $yllcorner + ($nrows-1-$i)*10;
    for ($j=0; $j<$ncols; $j++) {
      if (!$coarse_tile_mask[$i][$j]) { next; }
      # could also just set $v equal to $i...
      $v_by_lat{$lat} = int((80-$lat)/10);
      $lon = $xllcorner + 10*$j;
      $latmin = $lat;
      $latmax = $lat + 10;
      $lonmin = $lon;
      $lonmax = $lon + 10;
      ($hmin, $hmax, $vmin, $vmax) = &compute_hv_bounds_around_latlon_box($latmin,$latmax,$lonmin,$lonmax,$pix_per_deg_y);
      $hmin_by_latlon{$lat}{$lon} = $hmin;
      $hmax_by_latlon{$lat}{$lon} = $hmax;
      @tmplist = ();
      for ($h=$hmin; $h<=$hmax; $h++) {
        push @tmplist, $h;
      }
      $h_by_latlon{$lat}{$lon} = @tmplist;
      foreach $h (@tmplist) {
        $valid_hs_by_v{$v_by_lat{$lat}}{$h} = 1;
      }
    }
  }
}
else {
  # define hvs for the specified latlon box
  ($hmin, $hmax, $vmin, $vmax) = &compute_hv_bounds_around_latlon_box($latmin_in,$latmax_in,$lonmin_in,$lonmax_in,$pix_per_deg_y);
  $v_by_lat{$latmin_in} = $vmin;
  $hmin_by_latlon{$latmin_in}{$lonmin_in} = $hmin;
  $hmax_by_latlon{$latmin_in}{$lonmin_in} = $hmax;
  @tmplist = ();
  for ($h=$hmin; $h<=$hmax; $h++) {
    push @tmplist, $h;
  }
  $h_by_latlon{$latmin_in}{$lonmin_in} = @tmplist;
  foreach $h (@tmplist) {
    $valid_hs_by_v{$v_by_lat{$latmin_in}}{$h} = 1;
  }
}
  
# loop over 10-deg bands (modis v) from South to North
for ($i=$nrows-1; $i>=0; $i--) {
  if ($coarse_tile_mask_file ne "null") {
    $lat = $yllcorner + ($nrows-1-$i)*10;
  }
  else {
    $lat = $latmin_in;
  }
  $v = $v_by_lat{$lat};
  $vstr = sprintf "%02d", $v;

  if ($lat < $latmin_in || $lat > $latmax_in) { next; }

  # loop over 10x10 tiles
  for ($j=0; $j<$ncols; $j++) {
    if ($coarse_tile_mask_file ne "null") {
      if (!$coarse_tile_mask[$i][$j]) { next; }
      $lon = $xllcorner + 10*$j;
    }
    else {
      $lon = $lonmin_in;
    }

    if ($lon < $lonmin_in || $lon > $lonmax_in) { next; }

    # download MODIS files and aggregate
    if ($coarse_tile_mask_file ne "null") {
      $latmax = $lat+10;
      $lonmax = $lon+10;
    }
    else {
      $latmax = $latmax_in;
      $lonmax = $lonmax_in;
    }
    $script = "wrap_download_join_and_agg_MODIS_over_landcover.pl";
    $logfile = "log.$script.$domain.$lcid.$lat.$latmax.$lon.$lonmax.txt";
    $stagelist = "1,2";
    $cmd = "$script $lcformat $lcdir $lcpfx $lcid $lc_table $startyear $endyear $v $v $hmin_by_latlon{$lat}{$lon} $hmax_by_latlon{$lat}{$lon} $lat $latmax $lon $lonmax $resolution_out $landmask $domain $outdir $outpfx $force $stagelist > $logfile 2> $logfile &";
    print "$cmd\n";
    (system($cmd)==0) or die "$0: ERROR: $cmd failed\n";
  }

}

sub compute_hv_bounds_around_latlon_box {

  my $latmin = shift @_;
  my $latmax = shift @_;
  my $lonmin = shift @_;
  my $lonmax = shift @_;
  my $pix_per_deg_y = shift @_;

  my $lat, $lon, $htmp, $vtmp, $row, $col;
  my $hmin = -1;
  my $hmax = -1;
  my $vmin = -1;
  my $vmax = -1;
  my $resolution = 1/$pix_per_deg_y;

  $latmin += 0.5*$resolution;
  $latmax -= 0.5*$resolution;
  $lonmin += 0.5*$resolution;
  $lonmax -= 0.5*$resolution;

  foreach $lat ($latmin,$latmax) {
    foreach $lon ($lonmin,$lonmax) {
      ($htmp,$vtmp,$row,$col) = &latlon_to_modis_sinusoidal($lat,$lon,$pix_per_deg_y);
      if ($hmin == -1) {
        $hmin = $hmax = $htmp;
        $vmin = $vmax = $vtmp;
      }
      else {
        if ($htmp < $hmin) { $hmin = $htmp; }
        if ($htmp > $hmax) { $hmax = $htmp; }
        if ($vtmp < $vmin) { $vmin = $vtmp; }
        if ($vtmp > $vmax) { $vmax = $vtmp; }
      }
    }
  }

  return ($hmin, $hmax, $vmin, $vmax);
}
